{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/Qiskit/qiskit-tutorials/master/images/qiskit-heading.png\" alt=\"Note: In order for images to show up in this jupyter notebook you need to select File => Trusted Notebook\" width=\"500 px\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum Tomography Overview\n",
    "\n",
    "***\n",
    "### Contributors\n",
    "Gadi Aleksandrowicz (gadia@il.ibm.com), Christopher J. Wood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import itertools\n",
    "import qiskit \n",
    "from qiskit import QuantumRegister, QuantumCircuit\n",
    "from qiskit import Aer\n",
    "import qiskit.ignis.verification.tomography as tomo\n",
    "from qiskit.quantum_info import state_fidelity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The General Theory\n",
    "\n",
    "Quantum tomography is an experimental procedure to reconstruct a description of part of quantum system from the measurement outcomes of a specific set of experiments. In Qiskit Ignis we are currently concerned with the following two tomography tasks:\n",
    "\n",
    "1. **Quantum state tomography**: Given a state-preparation circuit that prepares a system in a state, reconstruct a description of the density matrix $\\rho$.\n",
    "2. **Quantum process tomograhpy**: Given a circuit, reconstruct a description of the quantum channel $\\mathcal{E}$ that describes the circuit.\n",
    "\n",
    "In both cases we rely on the assumption that we have access to a large number of identical copies of the system and so can perform several different measures on it.\n",
    "\n",
    "We can roughly split the tomography process to three stages:\n",
    "1. Preperation: Add suitable initialization/measurement devices to the quantum system.\n",
    "2. Experiment: Obtain measurement data from the quantum system.\n",
    "3. Tomography: Use the obtained data to reconstruct the system's description.\n",
    "\n",
    "Steps 1 and 2 are related to the quantum system being studied, whereas step 3 is a classical computation which can be carried out on standard computers.\n",
    "\n",
    "## State Tomography Overview\n",
    "\n",
    "Quantum state tomography is a method of reconstructing a description of the quantum state of a system from a set of experiments. While the state of ideal quantum system is described by a state-vector the state of an open quantum system (one that may experiences noise or other errors) is given by a density matrix $\\rho$. Quantum state tomography aims to reconstruct this density matrix. \n",
    "\n",
    "To do this we assume that the state $\\rho$ can be reliably prepared by a state-preparation circuit, and that itcan be subjected to several measurements with respect to different operators; this data can be used to reconstruct $\\rho$ or a close approximation of it by several different methods.\n",
    "\n",
    "\n",
    "### Definitions\n",
    "\n",
    "We denote by $\\mathcal{X}$ the state space of a closed (ideal) quantum system. In quantum computing this is typically the tensor product of $N$ 2-dimensional (qubit) systems $\\mathcal{X} = \\mathbb{C^{2N}}$. Valid quantum states $|\\psi\\rangle \\in \\mathcal{X}$ are those with norm-1:$|\\langle\\psi|\\psi\\rangle|^2 = 1$. \n",
    "\n",
    "We denote by $L(\\mathcal{X})$ the state space of linear maps on $\\mathcal{X}$, ($L: \\mathcal{X}\\rightarrow\\mathcal{X}$). The density-matrix for quantum system with state space $\\mathcal{X}$ is a linear map $\\rho \\in L(\\mathcal{X})$ that is also positive-semidefinite, and has trace equal to 1:\n",
    "1. **Unit trace:** $\\text{tr}[\\rho] = 1$\n",
    "2. **Positive-semidefinite:**: For all $|\\psi\\rangle \\in \\mathcal{X}$, $\\langle\\psi|\\rho|\\psi\\rangle \\ge 0$. This is denoted by $\\rho \\ge 0$.\n",
    "\n",
    "\n",
    "\n",
    "### Example: 1-qubit reconstruction using the Pauli basis\n",
    "\n",
    "Given the Pauli matrices $\n",
    "I=\\left(\\begin{array}{cc}\n",
    "1 & 0\\\\\n",
    "0 & 1\n",
    "\\end{array}\\right),\n",
    "X=\\left(\\begin{array}{cc}\n",
    "0 & 1\\\\\n",
    "1 & 0\n",
    "\\end{array}\\right),\n",
    "Y=\\left(\\begin{array}{cc}\n",
    "0 & -i\\\\\n",
    "i & 0\n",
    "\\end{array}\\right),\n",
    "Z=\\left(\\begin{array}{cc}\n",
    "1 & 0\\\\\n",
    "0 & -1\n",
    "\\end{array}\\right)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = np.array([[1,0],[0,1]])\n",
    "X = np.array([[0,1],[1,0]])\n",
    "Y = np.array([[0,-1j],[1j,0]])\n",
    "Z = np.array([[1,0],[0,-1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is easy to see they constitute an orthonormal basis for $M_2(\\mathbb{C})$ with respect to the Hilbert-Schmidt inner product $\\left\\langle A,B\\right\\rangle =\\frac{1}{2}\\text{tr}\\left(B^{\\dagger}A\\right)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def HS_product(A,B):\n",
    "    return 0.5*np.trace(np.conj(B).T @ A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "And hence,\n",
    "\n",
    "$$ \\rho =\\left\\langle \\rho,I\\right\\rangle I+\\left\\langle \\rho,X\\right\\rangle X+\\left\\langle \\rho,Y\\right\\rangle Y+\\left\\langle \\rho,Z\\right\\rangle Z = $$\n",
    "\n",
    "$$=\\frac{\\text{tr}\\left(\\rho\\right)I+\\text{tr}\\left(X\\rho\\right)X+\\text{tr}\\left(Y\\rho\\right)Y+\\text{tr}\\left(Z\\rho\\right)Z}{2}$$\n",
    "\n",
    "The values of $\\text{tr}\\left(X\\rho\\right), \\text{tr}\\left(Y\\rho\\right), \\text{tr}\\left(Z\\rho\\right)$ are the expectation values of $X$, $Y$, $Z$, respectively, and can be approximated by repeated measuring in the $X, Y$ and $Z$ bases. Since $\\text{tr}\\left(\\rho\\right)=1$ there is no need for additional measurements for the coefficient of $I$.\n",
    "\n",
    "### Example: 1-qubit Linear inversion\n",
    "\n",
    "The above method can be rephrased in more general form. First, any hermitian operator $H$ has a spectral decomposition of the form $H=\\sum \\lambda_i P_i$ where $\\lambda_i$ is an eigenvalue of $H$ and $P_i$ is the projection operator to the corresponding eigenspace. For the hermitian operators $X,Y,Z$ whose eigenvalues are 1 and -1 we can therefore write\n",
    "\n",
    "* $X = P^X_0-P^X_1$\n",
    "* $Y = P^Y_0-P^Y_1$\n",
    "* $Z = P^Z_0-P^Z_1$\n",
    "\n",
    "Where\n",
    "\n",
    "\n",
    "\n",
    "$$P^X_0=\\frac{1}{2}\\left(\\begin{array}{cc}\n",
    "1 & 1\\\\\n",
    "1 & 1\n",
    "\\end{array}\\right), P^X_1=\\frac{1}{2}\\left(\\begin{array}{cc}\n",
    "1 & -1\\\\\n",
    "-1 & 1\n",
    "\\end{array}\\right)$$\n",
    "\n",
    "$$P^Y_0=\\frac{1}{2}\\left(\\begin{array}{cc}\n",
    "1 & -i\\\\\n",
    "i & 1\n",
    "\\end{array}\\right), P^Y_1=\\frac{1}{2}\\left(\\begin{array}{cc}\n",
    "1 & i\\\\\n",
    "-i & 1\n",
    "\\end{array}\\right)$$\n",
    "\n",
    "$$P^Z_0=\\left(\\begin{array}{cc}\n",
    "1 & 0\\\\\n",
    "0 & 0\n",
    "\\end{array}\\right), P^Z_1=\\left(\\begin{array}{cc}\n",
    "0 & 0\\\\\n",
    "0 & 1\n",
    "\\end{array}\\right)$$\n",
    "\n",
    "In the Ignis code, these matrices are defined in **tomography.fitters.utils.pauli_preparation_matrix**. We give an explicit definition here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "PX0 = 0.5*np.array([[1, 1], [1, 1]])\n",
    "PX1 = 0.5*np.array([[1, -1], [-1, 1]])\n",
    "\n",
    "PY0 = 0.5*np.array([[1, -1j], [1j, 1]])\n",
    "PY1 = 0.5*np.array([[1, 1j], [-1j, 1]])\n",
    "\n",
    "PZ0 = np.array([[1, 0], [0, 0]])\n",
    "PZ1 = np.array([[0, 0], [0, 1]])\n",
    "\n",
    "projectors = [PX0, PX1, PY0, PY1, PZ0, PZ1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By Born's rule, $\\text{tr}\\left(P_{i}^{X}\\rho\\right)$ is the probability for the outcome $\\left|i\\right\\rangle$ when measuring in the X-basis, and this probability can be estimated directly using repeated meausrements in the X-basis. The $Y$ and $Z$ bases are handled\n",
    "similarily.\n",
    "\n",
    "The computation $\\text{tr}\\left(P_{i}^{X}\\rho\\right)$ can be replaced by the scalar product $\\vec{P}_i^X \\cdot \\vec{\\rho}$ where $\\vec{E}$ denotes the vector obtained from the operator $E$ by flattening its matrix (the result vector consists of the first row, then the second row etc.)\n",
    "\n",
    "Now we can construct a matrix $$M=\\left(\\begin{array}{c}\n",
    "\\vec{P}_{0}^{X}\\\\\n",
    "\\vec{P}_{1}^{X}\\\\\n",
    "\\vec{P}_{0}^{Y}\\\\\n",
    "\\vec{P}_{1}^{Y}\\\\\n",
    "\\vec{P}_{0}^{Z}\\\\\n",
    "\\vec{P}_{1}^{Z}\n",
    "\\end{array}\\right)$$\n",
    "\n",
    "Such that $$M\\vec{\\rho}=\\vec{p}=\\left(\\begin{array}{c}\n",
    "p_{\\left|0\\right\\rangle }^{X}\\\\\n",
    "p_{\\left|1\\right\\rangle }^{X}\\\\\n",
    "p_{\\left|0\\right\\rangle }^{Y}\\\\\n",
    "p_{\\left|1\\right\\rangle }^{Y}\\\\\n",
    "p_{\\left|0\\right\\rangle }^{Z}\\\\\n",
    "p_{\\left|1\\right\\rangle }^{Z}\n",
    "\\end{array}\\right)$$\n",
    "\n",
    "Is the equation relating the density operator to the observed probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.5+0.j ,  0.5+0.j ,  0.5+0.j ,  0.5+0.j ],\n",
       "       [ 0.5+0.j , -0.5+0.j , -0.5+0.j ,  0.5+0.j ],\n",
       "       [ 0.5+0.j ,  0. -0.5j,  0. +0.5j,  0.5+0.j ],\n",
       "       [ 0.5+0.j ,  0. +0.5j,  0. -0.5j,  0.5+0.j ],\n",
       "       [ 1. +0.j ,  0. +0.j ,  0. +0.j ,  0. +0.j ],\n",
       "       [ 0. +0.j ,  0. +0.j ,  0. +0.j ,  1. +0.j ]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = np.array([p.flatten() for p in projectors])\n",
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $M$ can be computed by knowing the operators used in the tomography, and the vector $\\vec{p}$ of probabilities can be estimated using the tomography results, all that remains is to solve the equation $M\\vec{\\rho}=\\vec{p}$ for $\\vec{\\rho}$. If the rank of $M$ is large enough this can be done by multiplying both sides by $M^\\dagger$:\n",
    "\n",
    "$M^\\dagger M\\vec{\\rho} = M^\\dagger \\vec{p}$\n",
    "\n",
    "$\\vec{\\rho} = (M^\\dagger M)^{-1} M^\\dagger \\vec{p}$\n",
    "\n",
    "In our example, we obtain the matrix \n",
    "$$(M^\\dagger M)^{-1} M^\\dagger = \\left(\\begin{array}{cccccc}\n",
    "\\frac{1}{6} & \\frac{1}{6} & \\frac{1}{6} & \\frac{1}{6} & \\frac{4}{6} & -\\frac{2}{6}\\\\\n",
    "\\frac{1}{2} & -\\frac{1}{2} & \\frac{i}{2} & -\\frac{i}{2} & 0 & 0\\\\\n",
    "\\frac{1}{2} & -\\frac{1}{2} & -\\frac{i}{2} & \\frac{i}{2} & 0 & 0\\\\\n",
    "\\frac{1}{6} & \\frac{1}{6} & \\frac{1}{6} & \\frac{1}{6} & -\\frac{2}{6} & \\frac{4}{6}\n",
    "\\end{array}\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "M_dg = np.conj(M).T\n",
    "linear_inversion_matrix = np.linalg.inv(M_dg @ M) @ M_dg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplication by the linear inversion matrix performs the reconstruction stage described earlier to obtain the density operator.\n",
    "\n",
    "### Example: 2-qubit Linear inversion\n",
    "\n",
    "For multiple qubit systems the technique of linear inversion remains the same. The projector operators are tensor products of 1-qubit projectors: $6^n$ projectors in total, since we measure according to $3^n$ operators (tensor products of $X,Y,Z$) and each operator has two projectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "projectors_2 = [np.kron(p1, p2) for (p1, p2) in itertools.product(projectors, repeat = 2)]\n",
    "M_2 = np.array([p.flatten() for p in projectors_2])\n",
    "M_dg_2 = np.conj(M_2).T\n",
    "linear_inversion_matrix_2 = np.linalg.inv(M_dg_2 @ M_2) @ M_dg_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now attempt to reconstruct the Bell state $\\frac{\\left|00\\right\\rangle +\\left|11\\right\\rangle }{\\sqrt{2}}$ from simulated tomography results. First, we prepare a quantum circuit which generates this bell state from the input $\\left|00\\right\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'OPENQASM 2.0;\\ninclude \"qelib1.inc\";\\nqreg q0[2];\\nh q0[0];\\ncx q0[0],q0[1];\\n'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q2 = QuantumRegister(2)\n",
    "bell = QuantumCircuit(q2)\n",
    "bell.h(q2[0])\n",
    "bell.cx(q2[0], q2[1])\n",
    "bell.qasm()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use Ignis' **state_tomography_circuits** procedure which generates the $3^n$ circuits obtained by adding to the bell circuit a measurement according to each of our measurement operators (Pauli by default). Then we execute on a standard simulator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "qst_bell = tomo.state_tomography_circuits(bell, q2)\n",
    "job = qiskit.execute(qst_bell, Aer.get_backend('qasm_simulator'), shots=5000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we load the data into the **StateTomographyFitter** which takes results data and can fit to a density matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "statefit = tomo.StateTomographyFitter(job.result(), qst_bell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the data we loaded into the **StateTomographyFitter** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{('X', 'X'): {'00': 2492, '11': 2508},\n",
       " ('X', 'Y'): {'00': 1256, '11': 1245, '10': 1262, '01': 1237},\n",
       " ('X', 'Z'): {'00': 1283, '11': 1286, '10': 1161, '01': 1270},\n",
       " ('Y', 'X'): {'00': 1253, '11': 1277, '10': 1263, '01': 1207},\n",
       " ('Y', 'Y'): {'10': 2478, '01': 2522},\n",
       " ('Y', 'Z'): {'00': 1311, '11': 1283, '10': 1198, '01': 1208},\n",
       " ('Z', 'X'): {'00': 1261, '11': 1207, '10': 1256, '01': 1276},\n",
       " ('Z', 'Y'): {'00': 1243, '11': 1251, '10': 1262, '01': 1244},\n",
       " ('Z', 'Z'): {'00': 2500, '11': 2500}}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "statefit.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use a private function **\\_fitter_data** to explicitly extract the probability vector $\\vec{p}$ and projector matrix $M$ that satisfy $M\\vec{\\rho} = \\vec{p}$. For typical usage we don't need to expose this data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "p, M, weights = statefit._fitter_data(True, 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use the linear inversion technique to reconstructo $\\vec{\\rho}$. Since we usually represent density matrices as matrices and not vectors, we use Numpy's **reshape** function to convert $\\vec{\\rho}$ into $\\rho$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 5.03133333e-01+0.j          5.36666667e-03+0.0095j\n",
      "  -3.56666667e-03-0.00053333j  5.00000000e-01+0.0031j    ]\n",
      " [ 5.36666667e-03-0.0095j      1.66666667e-03+0.j\n",
      "   2.77555756e-17-0.0029j      2.83333333e-03+0.00066667j]\n",
      " [-3.56666667e-03+0.00053333j  2.77555756e-17+0.0029j\n",
      "  -1.66666667e-03+0.j         -8.43333333e-03-0.0093j    ]\n",
      " [ 5.00000000e-01-0.0031j      2.83333333e-03-0.00066667j\n",
      "  -8.43333333e-03+0.0093j      4.96866667e-01+0.j        ]]\n"
     ]
    }
   ],
   "source": [
    "M_dg = np.conj(M).T\n",
    "linear_inversion_matrix = np.linalg.inv(M_dg @ M) @ M_dg\n",
    "rho_bell = linear_inversion_matrix @ p\n",
    "rho_bell = np.reshape(rho_bell, (4, 4))\n",
    "print(rho_bell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check the quality of our solution, we compute the fidelity between the real quantum state (obtained via simulation by a simulator that can return state vectors) and our calculated $\\rho$. The closer the fidelity to 1, the better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit Fidelity linear inversion = 0.9999999999999998\n"
     ]
    }
   ],
   "source": [
    "job = qiskit.execute(bell, Aer.get_backend('statevector_simulator'))\n",
    "psi_bell = job.result().get_statevector(bell)\n",
    "F_bell = state_fidelity(psi_bell, rho_bell)\n",
    "print('Fit Fidelity linear inversion =', F_bell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum Likelihood \n",
    "\n",
    "Linear inversion works perfectly on accurate data, but tomography data is never fully accurate. Two obvious obstacles are\n",
    "1. Since the number of measurements is limited, we do not obtain the probability vector $\\vec{p}$ but an approximation.\n",
    "2. The measurement process might be noisy.\n",
    "\n",
    "This may result in non-accurate or even self-contradicting $\\vec{p}$, and the result of linear inversion might not be a density function at all (e.g. not nonnegative, or trace different than 1).\n",
    "\n",
    "Since we want to solve the linear problem $A\\vec{x}=\\vec{p}$ for $x$, we can turn it into an optimization problem by attempting to minimize $\\|A\\vec{x}-\\vec{p}\\|_2$ while subjecting $x$ to additional constraints to ensure it is indeed a density matrix. This is done by **state_cvx_fit**.\n",
    "\n",
    "Another approach is to solve this optimization problem with no further constraints. The result might not be a density operator, i.e. positive semidefinite with trace 1; in this case the algorithm first rescales in order to obtain a density operator. This is done using **state_mle_fit**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit Fidelity CVX fit = 0.9998549234477454\n",
      "Fit Fidelity MLE fit = 0.9939629890765465\n"
     ]
    }
   ],
   "source": [
    "rho_cvx_bell = statefit.fit(method='cvx')\n",
    "F_bell = state_fidelity(psi_bell, rho_cvx_bell)\n",
    "print('Fit Fidelity CVX fit =', F_bell)\n",
    "\n",
    "rho_mle_bell = statefit.fit(method='lstsq')\n",
    "F_bell = state_fidelity(psi_bell, rho_mle_bell)\n",
    "print('Fit Fidelity MLE fit =', F_bell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Process Tomography Overview\n",
    "\n",
    "A quantum process, or quantum channel, describes the fixed-time evolution of a density matrix $\\rho^\\prime = \\mathcal{E}(\\rho)$. While ideal quantum processes are described by unitary matrix evolution, general quantum process can include the effect of errors in the ideal evolution such as those due to unwanted interaction with an environment system, systematic errors in gate operations, or other error processes. \n",
    "\n",
    "Quantum process tomography aims to find a description of $\\mathcal{E}$ from a set of experiments. Like state tomography, we measure with respect to several different bases. Process tomography also uses several different initial states, as opposed to state tomography where only $\\left|0^n\\right\\rangle$ is used. It can be thought of as performing state-tomography on tomographically complete set of a complete set of input states for the system.\n",
    "\n",
    "\n",
    "### Definitions\n",
    "\n",
    "Let $D(\\mathcal{X})$ be the set of density-matrices on a quantum system. Quantum channels are maps $\\mathcal{E}:L(\\mathcal{X})\\rightarrow L(\\mathcal{X})$ that are *completely positive-trace preserving* (CPTP):\n",
    "\n",
    "1. **CP:** $\\mathcal{E}$ is *completely-positive* if and only if $(\\mathcal{I}\\otimes\\mathcal{E})(\\rho)\\ge 0$ for all $\\rho \\ge 0$, where $\\mathcal{I}$ is an identity channel on an ancilla system.\n",
    "2. **TP:** $\\mathcal{E}$ is *trace-preserving* if and only if $\\text{Tr}[\\mathcal{E}(\\rho)] = \\text{Tr}[\\rho]$ for all $\\rho$.\n",
    "\n",
    "Together these two requirements ensure that the output of the quantum process will always be a valid density matrix.\n",
    "\n",
    "There are numerous representations for a quantum channel $\\mathcal{E}$. The simplest way to perform process tomogrpahy is to reconstruct a description of $\\mathcal{E}$ known as the *Choi-matrix* or process-matrix. This is because it is analogous to a density matrix $\\rho$, and hence the reconstruction is similar to quantum state tomography. Another common representation is the *Kraus* representation which we will also describe and show how it is related to the Choi-matrix. Another common representation is the Pauli Tranfer Matrix (PTM), which is useful for evolution since compositing channels becomes matrix multiplication. For more details about different representations of quantum channels and their relationships see [this paper](https://arxiv.org/abs/1111.6950).\n",
    "\n",
    "#### Choi-matrix representation\n",
    "\n",
    "The Choi-matrix description of a quantum channel $\\mathcal{E}$ is a unique bipartite matrix $\\Lambda_{\\mathcal{E}}\\in L(\\mathcal{X}\\otimes\\mathcal{X})$ given by the *Choi-Jamiolkowski isomorphism*.\n",
    "Given a basis $\\left\\{\\left|0\\right\\rangle, \\left|1\\right\\rangle,\\dots, \\left|n-1\\right\\rangle\\right\\}$ of $\\mathcal{X}$, the Choi-matrix $\\Lambda_{\\mathcal{E}}$ is defined as:\n",
    "\n",
    "\n",
    "$$\\Lambda_\\mathcal{E} = \\sum_{i,j=0}^{n-1} \\left|i\\right\\rangle\\left\\langle j\\right| \\otimes \\mathcal{E}(\\left|i\\right\\rangle\\left\\langle j\\right|)$$\n",
    "\n",
    "Evolution of a quantum state can be described in terms of this matrix by \n",
    "$$\\mathcal{E}(\\rho) = \\text{Tr}_1[\\Lambda_{\\mathcal{E}}(\\rho^T\\otimes\\mathbb{1})]$$\n",
    "\n",
    "where $\\text{Tr}_1$ denotes **partial trace**: $\\text{Tr}_1[A\\otimes B] = \\text{Tr}[A]\\,B$ and $\\mathbb{1}$ is the identity operator on $X$.\n",
    "\n",
    "In terms of the Choi-matrix, the requirements that $\\mathcal{E}$ is CPTP are:\n",
    "\n",
    "1. **CP:** $\\mathcal{E}$ is *completely-positive* if and only if the Choi-matrix is positive-semidefinite ( $\\Lambda_{\\mathcal{E}} \\ge 0$).\n",
    "2. **TP:** $\\mathcal{E}$ is *trace-preserving* if and only if $\\text{Tr}_2[\\Lambda_{\\mathcal{E}}] = \\mathbb{id}$.\n",
    "\n",
    "where $\\text{Tr}_2$ denotes **partial trace**: $\\text{Tr}_2[A\\otimes B] =\\text{Tr}[B]\\,A$.\n",
    "\n",
    "#### Kraus representation\n",
    "\n",
    "Another commonly used description of a quantum channel $\\mathcal{E}$ is via **Kraus operators**. A set $\\left\\{K_1, K_2, \\dots, K_t\\right\\}$ such that $K_i \\in L(X)$. In terms of this representation the evolution of a state is given by \n",
    "\n",
    "$$\\mathcal{E}(\\rho) = \\sum_{i=1}^k K_i\\rho K_i^\\dagger$$\n",
    "\n",
    "In terms of the Kraus operators the requirement that $\\mathcal{E}$ is CPTP is equivalent to $\\sum_{i=1}^k K_i^\\dagger K_i = \\mathbb{1}$.\n",
    "\n",
    "This is a common description of quantum processes. In case of a unitary operator $U$, the set $\\left\\{U\\right\\}$ is a Kraus operator representation of it,\n",
    "$$\\mathcal{E}(\\rho) = U \\rho U^{\\dagger}$$\n",
    "Measurements and noises are also commonly given via Kraus operators. Note however that the Kraus representation of a quantum process is not unique.\n",
    "\n",
    "The Choi-matrix may be constructed from Kraus operators via \n",
    "$$\\Lambda_\\mathcal{E} = \\sum_{i,j=0}^{n-1}\\sum_{l=1}^k \\left|i\\right\\rangle\\left\\langle j\\right| \\otimes K_l\\left|i\\right\\rangle\\left\\langle j\\right|K_l^\\dagger = \\sum_{l=1}^k |K_l\\rangle\\!\\rangle\\!\\langle\\!\\langle K_l|$$\n",
    "\n",
    "where $|A\\rangle\\!\\rangle$ denotes a vectorized matrix, which is a column-vector obtained by stacking the columns of $A$.\n",
    "\n",
    "To construct a set of Kraus operators from a Choi-matrix we may use the spectral-decomposition of $\\Lambda_{\\mathcal{E}}$ and define $ |K_l\\rangle\\!\\rangle\\ = \\sqrt{\\lambda_l}|v_l\\rangle$, where ${\\lambda_l} \\ge 0$ and $|v_l\\rangle \\in \\mathcal{X}\\otimes\\mathcal{X}$ are the eigenvalues and eigenvectors respectively of $\\Lambda_{\\mathcal{E}} \\ge 0$.\n",
    "\n",
    "#### Pauli-Transfer-Matrix (PTM) representation\n",
    "\n",
    "Another description of a quantum channel $\\mathcal{E}$ is via **PTM**. The elements of the PTM are given as \n",
    "\n",
    "$$R_{ij} = \\frac{1}{d}\\text{Tr}[\\sigma_i\\mathcal{E}(\\sigma_j)] $$\n",
    "\n",
    "where $\\sigma$ are the Pauli matrices ($4^n$ possible matrices). The advantage of this method is that composing channels,\n",
    "\n",
    "$$\\mathcal{E}(\\rho) = \\mathcal{E_1}(\\mathcal{E_2}(\\rho)) $$\n",
    "\n",
    "is\n",
    "\n",
    "$$R = R_1 \\cdot R_2 $$\n",
    "\n",
    "### Process tomography with the Choi matrix\n",
    "\n",
    "The CP condition on the Choi-matrix ($\\Lambda_{\\mathcal{E}})\\ge 0$ is equivalent to the state-tomography condition that a density matrix is positive-semidefinite. The only difference is that the normalization of the Choi-matrix is not trace-1, but trace given by the dimension of $\\mathcal{X}$ ($\\text{Tr}[\\Lambda_{\\mathcal{E}}] = d$). This allows us to use the same fitters from state-tomography to fit the Choi-matrix adjusting for this difference in normalization. *The key difference between state tomography is the TP condition which is an additional constraint we must add to the tomography fitters if we want to ensure the fitted channel is CPTP.*\n",
    "\n",
    "To see this relationship explicity consider preparing the system in an initial state $\\rho$ and measuring. The probability of measurement for a projector $\\Pi$ is given by\n",
    "\n",
    "$$p_{ij} = \\text{Tr}[\\Lambda_{\\mathcal{E}}(\\rho^T\\otimes\\Pi_j)]$$\n",
    "\n",
    "Using the above equation we can see that the act of measuring the outcome of $\\mathcal{E}$ on some initial $\\rho$ with some projector $\\Pi$ can be seen as the act of measuring $\\Lambda_\\mathcal{E}$ (when it is considered as a **state** in the space $L(\\mathcal{X}\\otimes\\mathcal{X})$) with a measurement operator $\\overline{\\rho}\\otimes \\Pi$, where we have used that $\\rho^T=\\overline{\\rho}$ since $\\rho$ is Hermitian.\n",
    "\n",
    "This gives rise to the following algorithm for process tomography:\n",
    "\n",
    "1. Obtain a set of initial states $\\left\\{\\rho_1, \\dots \\rho_{k}\\right\\}$ and a set of projectors $\\left\\{P_1, \\dots P_{t}\\right\\}$ such that the set of all projectos $\\Pi_{i,j} = \\overline{\\rho_i}\\otimes P_j$ is tomographically complete for $L(\\mathcal{X}\\otimes\\mathcal{X})$.\n",
    "2. Obtain measurements of $\\Lambda_\\mathcal{E}$ with $\\Pi_{i,j}$ by initializing a system to $\\rho_{i}$, applying $\\mathcal{E}$ (e.g. via a simulator or the quantum computer to check) and measuring via $P_j$.\n",
    "3. Pass the results and the description of $\\Pi_{i,j}$ to an algorithm for state tomography and obtain $\\Lambda_\\mathcal{E}$. *Note that to ensure the fitted state is also TP requires adding the additional constraint to the fitter (if supported) that $\\text{Tr}_2[\\Lambda_{\\mathcal{E}}]=\\mathbb{1}$.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
